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A new class of multilayered functional mate¬ 
rials has recently emerged in which the compo¬ 
nent atomic layers are held together by weak van 
der Waals forces that preserve the structural in¬ 
tegrity and physical properties of each layer [T]. 
An exemplar of such a structure is a transis¬ 
tor device in which relativistic Dirac Fermions 
can resonantly tunnel through a boron nitride 
barrier, a few atomic layers thick, sandwiched 
between two graphene electrodes. An applied 
magnetic field quantises graphene’s gapless con¬ 
duction and valence band states into discrete 
Landau levels, allowing us to resolve individ¬ 
ual inter-Landau level transitions and thereby 
demonstrate that the energy, momentum and 
chiral properties of the electrons are conserved 
in the tunnelling process. We also demonstrate 
that the change in the semiclassical cyclotron 
trajectories, following a tunnelling event, is a 
form of Klein tunnelling for inter-layer transi¬ 
tions. 

An electron moving through the hexagonal crystal 
structure of graphene is not only quasi-relativistic but 
also exhibits chirality [2], which means that its wave- 
function amplitude is intrinsically coupled to the di¬ 
rection of motion. This gives rise to the phenomenon 
of Klein tunnelling whereby an electron can pass with 
unity transmission through a potential barrier formed in 
the graphene layer [am. In principle, chirality should 
affect the electronic properties of graphene-based de¬ 
vices. To investigate this effect we focus on a van der 
Waals heterostructure in which Dirac fermions can res¬ 
onantly tunnel between two graphene electrodes sepa¬ 
rated by a hexagonal boron nitride tunnel barrier [Ms]. 
Recent work on this type of transistor has demonstrated 
that even a small misalignment of the crystalline lat¬ 
tices of the two graphene electrodes lowers the trans¬ 
lation symmetry in the plane of the tunnel barrier and 
gives rise to an impulse which modifies the dynamics of 
the tunnelling electron CHia. By applying a quantising 
magnetic field perpendicular to the layers, we show that 
electron tunnelling is governed by the laws of conserva¬ 
tion of energy and of in-plane momentum. In addition, 
we find that the effect of electron chirality on the tunnel 
current is enhanced by a quantising magnetic field. We 



Figure 1: a Schematic of the device showing the two mis¬ 
aligned graphene lattices (bottom, red and top, blue) sep¬ 
arated by a boron nitride tunnel barrier, yellow, b dashed 
black lines show the Brillouin zone boundary for electrons 
in the bottom graphene layer. Red arrows show the vector 
positions of the Dirac points Kb^ (red circles) relative to 
the r point. Blue arrows show the positions of the Dirac 
points in the top layer, Kt^ (blue circles), misorientated at 
an angle 0 to the bottom layer. 

also demonstrate that, following an electron tunnelling 
transition, the semiclassical cyclotron trajectory of the 
electron changes in a way that is analogous to intra¬ 
layer Klein tunnelling. 

Our device, with bias, 14, and gate, Vg, voltages ap¬ 
plied, is shown schematically in Eig It consists 
of a 4-layer thick hexagonal boron nitride (hBN) tun¬ 
nel barrier m sandwiched between two high purity 
crystalline graphene electrodes. The doped Si layer 
of a Si 02 /n-Si substrate acts as the gate electrode. 
The two graphene lattices are intentionally aligned to 
within an angle of 1°, see ref. |8] for details. How¬ 
ever, even this slight misalignment, or “twist angle”, 0, 
leads to a significant /c—space displacement of magni- 
tude, AK = |AK=^| = |K± - K±| = 2sin(6»/2)|K±| 
of the Dirac cones at the corners of the Brillouin zones 
[HHiT], see Eig. and Eig. [^, insets. This displace¬ 
ment induces an impulse on tunnelling electrons and 
has a large effect on the measured current-voltage char¬ 
acteristics and their magnetic field dependence. 

Eig. 1^ (black dashed curve) shows the measured 
current-volt age curve, /( 14 ), at I 4 = 0 in the absence 
of a magnetic field. The current increases at a thresh¬ 
old bias voltage Vi and reaches a resonant peak when 
Vb = V 2 = O.bS y, beyond which there is a region of 
negative differential conductance. When I 4 = hi, see 
inset i of Eig. |^, the Eermi circle in one cone partially 
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Figure 2: ql I{Vh) curves measured when = 0 for B — ^ 
(black dashed) and 4 T (red solid), the latter offset by 7.5 
fiX for clarity. Insets i and ii show the relative energies of 
the displaced Dirac cones in /c—space, in the bottom (red) 
and top (blue) electrodes, whose intersections are shown by 
the thick yellow curves, at the voltages Vi and V 2 marked 
by the labelled vertical arrows. The Fermi circles of the two 
layers are shown in white, b differential conductance, G(Vb), 
measured at I 4 = — 40 V (blue lower curve) and Vg = 40 V 
(red upper curve) when B = 4 T and temperature T = 4 K. 
Upper curve is offset by 250 juS i.e. dotted lines mark G = 0 
for the two curves. 


overlaps with empty states in the other, so that elec¬ 
trons can tunnel with energy and momentum conserva¬ 
tion [SI [9] . When T4 = V 2 (inset ii) the cones intersect 
along a straight line and the current reaches a resonant 
maximum. 


A magnetic field, B, applied perpendicular to the 
graphene layers quantises the electron energy into a 
spectrum of unequally-spaced Landau levels (LLs) de- 
fined by ^ = sgn(nb^t)^/2lni^hvF/lB, where rib^t 
is the LL index in the bottom (b) and top (t) layers, 
and IB = ^JhjeB [181 EDHSQI [S^ • By comparing our 
measured tunnel current with transfer Hamiltonian cal¬ 
culations, we demonstrate the composite spatial-spinor 
form of the quantised Landau states and the effect of 
chirality on the measured current-voltage characteris¬ 
tics. In addition, by using a semiclassical description of 
the cyclotron orbits of an electron before and after the 
tunnelling event, we explain how resonant tunnelling is 
enabled by the large momentum impulse induced by the 
small twist angle between the two graphene lattices. 


a B = 2J, experiment 


b B = 2J, theory 
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Figure 3: Colour maps showing G{Vb,Vg) at T = 4 K mea¬ 
sured (a) and calculated (b) when H = 2 T and (c measured, 
d calculated) when H = 4 T. Colour scales for a,c are in 
/iS and for b,d normalised to the maximum conductance in 
the maps. Black and white dashed curves enclose regions 
around 14 = 0 within which only conduction-conduction 
(upper region with Vg >0), or only valence-valence (lower 
region with Vg < 0) tunnelling occurs. 


Effect of a perpendicular magnetic field on 
resonant tunnelling: experiment and theory 


Landau level quantisation induces weak features in 
I{Vb) when = 0 for 0.08 V < I 4 < 0.35 V (see region 
of the red curve in Fig. [^ indicated by the green hor¬ 
izontal bar) and sharp, large amplitude, resonant fea¬ 
tures in the differential conductance, G{Vb) = dl/dVb, 
as shown in Fig. for gate voltages Vg = ±40 V. 
By combining similar plots at intermediate gate volt¬ 
ages, we generate the colour maps of G{Vb,Vg) shown 
in Figs. and c, for H = 2 and 4 T, respectively. The 
regions of high conductance are patterned by small “is¬ 
lands” that originate from resonant tunnelling of elec¬ 
trons when LLs in the two graphene layers become 
aligned in energy (shown schematically in Fig. 

These islands are sharply defined close to 44 = 0 but 
become broadened at high |I4|, which could arise from 
carrier heating due to high current levels and/or in¬ 
creased lifetime broadening. 

We model our data (see Fig. [^,d) using a Bardeen 
transfer-Hamiltonian approach, taking the full two com¬ 
ponent form of the LL eigenstates and the following 
device parameters: the doping densities in the bottom 
and top graphene layers are 2.0 x 10^^ cm“^ (p-fype) 
and 3.6 x 10^^ cm“^ (n-type) respectively, and the twist 
angle 0 = 1°. A fit to the /(I4, Vg) curves at H = 0 pro¬ 
vides accurate values of these parameters [8] (also see 
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Figure 4: a,b Dirac cones showing the energy-wavevector dispersion relation, ^(k), for electrons in the bottom (red) and 
top (blue) graphene layers when S = 0 and 14 = 0.28 V a and 0.58 V b. Rings of constant energy on the surface of the 
cones show the energies and semiclassical /c—space radii of LLs with indices nt and rit. The black rings in a and b highlight 
= 1 to nt = 3 and n?, = 2 to = 16 transitions, respectively. Occupied electron states in the bottom (top) layer are 
shaded dark red (blue) up to the Fermi level, in that layer, c Colour map showing tunnelling rates, kF(nb,nt), for 
scattering-assisted transitions (taking (7 = 9 nm) between LLs with indices rih and rit in the bottom and top electrodes. The 
dotted and solid curves show the loci calculated using Eq. it- For all panels, B — AT. 


[HIES] and Supplementary Information, SI, for further 
details). 

Our model gives a good fit to the magneto-tunnelling 
data, in particular the shape and relative strength of 
the islands of high conductance. It enables a detailed 
analysis of the pattern of conductance peaks (see SI). 
We now focus on the underlying physics that controls 
the overall pattern of peak amplitudes, in particular 
the effect of twist angle and chirality on the tunnelling 
process. 

Transition rates between chiral LL eigenstates 

The displacement, AiC, of the Dirac cones due to 
the twist angle is shown schematically in Figs. I^,b. 
It can be represented by, and is equivalent to, the ef¬ 
fect of a strong pseudo-magnetic field applied parallel 
to the graphene layers m- We describe the combined 
effects of the misalignment and the Landau-quantising 
applied magnetic field by a vector potential in the Lan¬ 
dau gauge, 

Ab,( = {ihAK^, -eBx + IhAK^ ,Q) /e, (1) 

where I = 0,1 for the b, t layers. In a perpendicular 
magnetic field, the electron wavefunctions at the 
point have the analytic forms [20l [21] 


an equivalent contribution to the tunnelling matrix ele¬ 
ment, see SI. The tunnel rates between LLs, LF(n 5 ,nt), 
depend on the overlap integrals of the initial and fi¬ 
nal wavefunctions summed over the /c—states in the two 
layers (see SI) and therefore permits tunnelling between 
SHO states with a range of different n indices. Fig. Ilh,b 
show the energies and semiclassical trajectories (yellow 
rings) of the quantised Landau states. 

Fig# is a colour map of the inter-LL transition rate 
W{ni),nt) at 5 = 4 T (see Eq. (25) of the SI). It reveals 
narrow yellow regions where W{niy, rit) is high. In other 
areas (black), tunnelling is suppressed. The regions of 
high W{ni),nt) originate from the spatial form and rel¬ 
ative positions of the wavefunctions in the bottom and 
top electrodes. Within the upper right and lower left 
quadrants of the colour map, transitions between equiv¬ 
alent bands (conduction-conduction, c-c, and valence- 
valence, v-v) are strongly enhanced compared to tun¬ 
nelling between different bands (c-v and v-c). This 
asymmetry, found for all values of 5, is a consequence 
of chirality. In contrast, when we remove the effect of 
chirality from our model by using pure (single compo¬ 
nent) LL wavefunctions, the tunnelling matrix elements 
are the same for transitions between equivalent and dif¬ 
ferent bands (see SI). 


Effect of chirality on tunnel current 


exp(i%) 


-sgn{nby)i(l)\ 



( 2 ) 


The two-component chiral states comprise plane waves 
along y and simple harmonic oscillator (SHO) waves, 
0, along X with indices that differ by 1. The centres 
of the SHO wavefunctions in the top and bottom lay¬ 
ers are shifted by 1‘^AK^ and there is an additional 
plane wave factor for the top layer whose argument is 
AK+{x - Xt), where = ll{k + AX+). The Bloch 
states near the K~ point have similar form and make 


The asymmetry in the transition rate colour map in 
Fig# manifests itself in the observed pattern of con¬ 
ductance peak amplitudes. In certain regions of the 
{yb^yg) plot, tunnelling is exclusively between equiv¬ 
alent bands, as shown in Fig. Here, the black 
and white dashed curves bound the regions of 14 — 
space where tunnelling is either only c-c (upper region, 
> 0) or v-v (lower region, < 0), respectively. 
Within these regions the amplitudes of the resonant 
peaks are high, i.e. dark red. Increasing 14 beyond 
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the lower region induces a changeover from tunnelling 
between equivalent bands to a mixture of tunnelling be¬ 
tween equivalent and different bands and is therefore ac¬ 
companied by a suppression of the conductance peaks. 
This is a direct manifestation of electron chirality. This 
changeover also occurs as Vt, decreases across the left 
hand edge of the upper bounded region of Fig. Eh-d. 

The effect of chirality on the peak amplitudes in these 
regions is seen more clearly in the enlarged lower re¬ 
gion of the G{Vb,Vg) maps at 5 = 2 T shown in Figs, 
[^-c. In both our experiment, a, and calculations, b, 
the conductance peak amplitudes are larger within the 
bounded region in the lower left-hand side of the plot, 
labelled L in Fig. ii, where v-v tunnelling dominates 
and smaller in the bounded region in the lower right- 
hand side of the plot where tunnelling is a mixture of 
v-v and v-c transitions (region labelled R in Fig. |^). 
For comparison, in Fig. Ef we show G{Vb^Vg) calculated 
when chirality is “switched off”, i.e. with each eigen¬ 
state represented by a single SHO wavefunction with no 
pseudospin component (see SI). In contrast to the chiral 
theory and experimental data, the conductance peaks 
for the non-chiral calculations have similar amplitudes 
in regions L (v-v) and R (v-v and v-c). 

To quantify the effect of chirality on the tunnel cur¬ 
rent, we calculate the ratio of the mean conductance in 
region L to that in region R, {G)l/{G)r (see Fig. |^). 
In the bar chart in Fig. we show {G)l/{G)r when 
B = 0, 2 and 4 T. For each field value, {G)l/{G)r for 
the measured data (red) and the chiral calculations (yel¬ 
low) are similar to each other. In contrast {G)l/{G)r 
is significantly smaller for the non-chiral calculations 
(blue). In addition, with increasing B the difference be¬ 
tween the chiral and non-chiral results becomes larger: 
at higher B there are fewer LL transitions within regions 
L and R and, for those transitions that do occur, the dif¬ 
ference between the chiral and non-chiral conductance 
is more pronounced. Hence, the measured dependence 
of the conductance peak amplitudes on Vg^ 14, and H, 
reveals and demonstrates the chiral nature of the elec¬ 
trons and the associated asymmetry in the tunnelling 
rates (see Fig. Ef)- 


Nested and figure of 8 cyclotron orbits 

A semiclassical picture, in which electrons undergo 
cyclotron motion in both real- and /c-space, provides 
further insights into the physics of tunnelling in these 
devices. In /c-space, the orbital radii 
in the two graphene layers are separated by AK^. The 
solid and dotted curves in Fig. are loci of initial and 
final states along which the corresponding semiclassical 
orbits just touch, so that the tunnelling electrons can 
make a continuous classical trajectory in the (kx^ky) 
plane. These loci are defined by 

(3) 




Figure 5: a-c Colour maps showing G(Vb,Vg) when R = 2 
T. a,b are enlargements of the lower parts of the colour maps 
in Fig. |7^,b respectively. Panel a shows experimental data 
(T = 4 K), b is calculated using the full model with chiral 
electrons, and c calculated using non-chiral wavefunctions 
i.e. comprising a single simple harmonic oscillator state. 
Colour bars in a and b,c are in jaS and normalised units, 
respectively. Solid curves in a-c enclose regions of the colour 
map where tunnelling is only v-v (labelled L in d) or a mix¬ 
ture of v-v and v-c (labelled R in d). Bar charts in e show 
the ratio, {G) l / {G)r, of the mean conductance in regions L 
and R (see d) for the measured data (red), and calculated 
for chiral (yellow) and non-chiral (blue) electrons. 


Here the — and + signs specify, respectively, cyclotron 
orbits that describe a “figure of 8” (F-8) and nested (N) 
form. Examples are shown by the projected circles in 
the lower parts of Figs. and b. The spatial varia¬ 
tion of the real (dark) and imaginary (light) components 
of the corresponding two-component LL wavefunctions 
are also shown {x axis re-scaled by 1/1% to enable com¬ 
parison between the k-space trajectories and the spa¬ 
tial form of the SHO wavefunctions). The maxima in 
the wavefunction amplitude are located at the turning 
points of the semiclassical orbit so that, when Eq. ® 
is satisfied, i.e. along the solid (dotted) locus in Eig.|^ 
for N (E-8) orbits, the wavefunction overlap integral is 
large. 

The N and E-8 semiclassical orbits determine the de¬ 
pendence of G on H, Vb and Vg. At the onset of current 
(see red curve and arrow labelled Vi in Eig.[^) the ener¬ 
getically aligned LLs correspond to semiclassical orbits 
with the E-8 form, see black rings in Eig. [^. Conse¬ 
quently the matrix elements are large, allowing tunnel 
current to flow. At the resonant current peak (see red 
curve and arrow labelled V 2 in Fig.|^) the Dirac cones 
just touch and their intersection is a straight line. As 
a result, all energetically aligned LLs have high matrix 
elements because all the corresponding semiclassical or¬ 
bits have either E-8 or N forms, see black and yellow 
rings in Eig. When Vb increases beyond the current 
peak, many LLs that become aligned energetically have 
cyclotron orbits that do not overlap spatially and so the 


Kt = AK ± Kb- 
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tunnelling matrix elements and current decreases. 

The semiclassical analysis also highlights the effect 
of the lattice misalignment on the electron dynamics. 
At the point of intersection of the = 1 and rit = 
3 orbits, the electron “back-scatters” in both /c-space 
and real space, making a 180° direction change where 
the orbits touch in the (kx-ky) plane, see lower part 
of Fig. This change in kinetic momentum at the 
intersection between the two orbits is induced by the 
impulse, arising from the misorientation of the 

two graphene layers and the associated vector potential, 
which acts like an in-plane pseudo-magnetic field for 
tunnelling electrons, see Eq. Q. 

As shown in Fig.|^, for the F-8 orbits, the tunnelling 
transition reverses the wavevector in the bottom and 
top electrodes, and kt, measured relative to the Dirac 
point of the two layers. In contrast, for N orbits the 
direction of the wavevector in the two electrodes is un¬ 
changed during tunnelling; only its magnitude changes 
(Fig.|^). 


Cyclotron orbits and Klein tunnelling 

In graphene, the chiral nature of an electron in the 
absence of a magnetic field can be expressed by the ex¬ 
pectation value of the pseudospin operator with respect 
to the eigenstate. For the valley this expectation 
value is {a) = 5(±cos(^, sin(/9), where is the polar di¬ 
rection of the wavevector. Therefore, in our semiclassi¬ 
cal model, for N orbits in both valleys (cr) is unchanged 
for equivalent band transitions but is rotated by 180° for 
transitions between different bands. In contrast, for F-8 
orbits (cr) is reversed for transitions between equivalent 
bands and unchanged for transitions between different 
bands. 

When (cr) is unchanged, the inter-layer tunnelling 
process bears an analogy with intra-layer Klein tun¬ 
nelling US]. The Klein paradox is predicted to oc¬ 
cur for electrons tunnelling through a barrier in planar 
graphene where unity transmission is expected when 
the pseudospin is conserved. In our device, the tun¬ 
nelling electron makes a “quantum jump” across the 
barrier; hence, the tunnelling rate can be high even if 
pseudospin is reversed, provided there is strong spatial 
overlap between the initial and final LL wavefunctions. 
However, as for the case of Klein tunnelling in planar 
graphene, the orientation of (cr) in the initial and final 
states determines the tunnelling rate. Physically this is 
due to the interference between the A and B sublattices 
of graphene (see Eq. (13) of the SI and [2]). In our ex¬ 
periments, resonant tunnelling is enabled by the twist 
of the graphene electrodes. This provides the impulse 
to induce the momentum and orbit centre change re¬ 
quired for energy- and k-conserving tunnel transitions 
with high matrix elements. In particular, our data in¬ 
dicate that the pseudospin of the electrons is conserved 
for the tunnelling transitions at the current peak. 



Figure 6: a,b Upper: vertical (horizontal) curves show 
the real (imaginary) parts of the real space electron wave- 
function in the bottom (red curves) and top (blue curves) 
graphene electrodes respectively with H = 4 T and a = I 
(red) and nt = 3 (blue) and b = 2 (red) and nt = 16 
(blue). The x axis is scaled by 1% for comparison with lower 
plots: circles show corresponding figure of 8 and nested cy¬ 
clotron orbits in k- space {kx,ky axes inset and direction of 
motion marked by arrows) with orbit centres separated by 
AK. The vertical black lines connecting upper and lower 
parts of the figure show the classical turning points. 

Conclusions 

We have investigated how LL quantisation of Dirac- 
Weyl Fermions reveals the effects of chirality on the res¬ 
onant tunnelling transitions in graphene-hBN-graphene 
heterostructures. Semiclassically, when the electron 
tunnelling trajectory takes the form of off-centred 
“nested” or “figure of 8” transitions, the pseudospin is 
either unchanged or undergoes a pseudospin-fiip tran¬ 
sition of 180°. At the resonant peak of our measured 
and calculated current-voltage curves the pseudospin is 
conserved for all transitions, in analogy with Klein tun¬ 
nelling in single-layer graphene. Analysis of the experi¬ 
mental data confirms that the Dirac-Weyl model for the 
electronic states of electrons in graphene provides an ac¬ 
curate description of the tunnel current flowing perpen¬ 
dicular to the plane of the barrier in these stacked van 
der Waals heterostructures, so-called “vertical” trans¬ 
port. Our results demonstrate that the chirality pro¬ 
vides an important contribution to the characteristics 
of graphene-based tunnelling devices, and should there¬ 
fore be taken into account when designing future elec¬ 
tronic components based on materials with Dirac-like 
energy spectra. 
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SUPPLEMENTARY INFORMATION 
MODEL 

The graphene lattices in our device are slightly misorientated by an angle 6 > 1° which results in a relative 

displacement in the positions of the Dirac points in K space, AK^ = {R{0) — 1 )K^, where R{0) is the rotation 
matrix. The label ± corresponds to the two inequivalent K points with positions given by = ±( 47 r/ 3 a, 0 ), 
where a = 2.46 Ais the lattice constant of graphene. The Dirac points in the bottom electrode are at and in 
the top electrode + AK^. The relative shift of the Dirac points is analogous to an in-plane magnetic field. 
Therefore, we describe the displacement of the K points using the following vector potential for electrons in the 
bottom and top layers, which also includes the effect of a magnetic field, B that is applied perpendicular to the 
graphene layers. 


Ab,t = {ihAK^ ,-eBx + IhAK^ , O) /e, (4) 

where I = 0,1 in the bottom (b) and top (t) layers respectively. The electron momentum takes the form p ^ p+eA, 
so that the effective mass Hamiltonian for Dirac electrons in graphene becomes 

TT^ _ n, f ^ ^ {Py 

^ V (Px + e^x,b,t) + ^ {Py + eAy,b,t) 0 

where vf = 10^ ms“^. The Hamiltonian has the form of a quantum harmonic oscillator so that the electron has 
discrete Landau energy levels given by 



= sgn{nb,t)\nb,t\^eBhvp, ( 6 ) 

where is an integer that labels the energy levels in the two electrodes, positive for electrons in the conduction 
band and negative in the valence band and 


1 (n > 0 ) 

sgn(n) = ^ 0 (n = 0 ) 

[-1 (n< 0 ). 

The electron wavefunctions at the two Dirac points are therefore 


( 7 ) 


and 




K+ 

'^b,t ^kb,t 


(r) 


Erib^t 

Vl 


exp {ikb,ty) 


-sgn{nb,t)i4>\ 

rib, 



where 




K- 

nb,t ■)kb,t 


(r) 


Vl 


exp {ikb,ty) 


sgn(nb,t)i(/)| 

ribjt 

4’\nb,t\ 



(8) 

( 9 ) 


where 


and 


Cn 


1 (n = 0 ) 

1/V2 (n^O) 


Vd = 


y'2\^i>\\nb\\y^lB 


exp 




rr\nb\ ’ 


^\nt\ ~ 




exp 


.^{x-Xtf-iAKt{x-Xt) 


F^^\[f-{x-Xb) 


( 10 ) 


( 11 ) 


( 12 ) 


Here Ib = ^/h/eB and is the nth order Hermite polynomial. The orbit centre in the bottom and top electrodes 
are given by X 5 = l%ky and Xt = respectively. The effect of the misorientation of the two 

graphene sheets is to shift the relative position of their orbit centres by /^AAT^ and introduce a phase difference 


of AK,{x-Xt). 
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Matrix element 


We assume that electrons can undergo elastic scattering which we describe using a Gaussian scattering potential: 


Vs{x,y) = Voe 


(13) 


where <j ~ 10 nm is the scattering length scale. The matrix element for tunnelling between the bottom and top 
electrodes is given by 


Mu= [ dV^;{r,z)VsMr,z). 
Jv 


(14) 


First we consider the integral in the ^ direction. We assume that the electron wavefunctions decay exponentially 
into the barrier regions so that the integral is a constant, equal to 


^ Kd 

“ D ■ 


(15) 


where d is the barrier width. We assume hz to be independent of energy to facilitate analysis of the current. For 
full analysis of different dependent models for see Ref. [5]. In the basis of Bloch wavefunctions [9] and [8], 
the matrix element is given by 


^bt {p^bi '^t'l^b'i^t) — kf) (1^51, \Tii \ ^ k}) ^ kf) ^ Ix (| I5 1 5 ^65 ) 

±isgn{nb)Ix{\nb\, \nt\ - 1, h, kt) + sgn(n6)sgn(nT)4(|n6| - 1, \nt\ - 1, k^, h)] 


(16) 


where Ix and ly are the overlap integrals of the wavefunctions along the x and y axes respectively. On first 
inspection, Eq. (16) appears to reveal that the matrix element is different for tunnelling between valleys 


(upper sign) compared to that between K~ valleys (lower sign). However, AK+ = — AK“ and, consequently, it 
can be shown that the matrix element for transitions between the same valleys are equivalent. Our matrix element 
does not explicitly include the cell-periodic parts of the Bloch functions, ?i"’^(r), where a and P label the two 
atoms in grahene’s unit cell. This is because for small relative rotations of the two layers, the spatial overlap 
integral of the cell-periodic parts of the wavefunction j {R{0 )y)u^^^ [r) are approximately equivalent for all 

combinations of a and and therefore will only have a small quantitative effect on the matrix element [9]. 


Overlap integrals for seattering assisted tunnelling 
The overlap integrals ly and Ix can be shown [32] to have following form: 


ly = \/^crexp (—A/c^cr^/2) , 

within which A/c = kij — kf. The overlap integral in the x direction, is given by: 


(17) 


{Rb-) — j-i Pbt {Rb-i ^b-) ^t) 

Mb 


\nh,nt\ 

E 

i=o 


f- 


nt ^ (l — 


(18) 


{2ay 



1 

1 

cr 

H 

'aT-lB {h + I\K+y 

-3 



1 

(M 

1 


(19) 


where a = 1 /(Ib, 
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within which 


^ = T2 - ^ ((fc* + AK^f + kl) - iAKt {h + AK^) , 

and 

T = ^^{h + AK^ + h + iAKt). 


( 22 ) 


(23) 


Current 

The current between the layers is given by the sum over states in the top and bottom layers: 

bt 

where the Fermi functions for the bottom and top layers are given, respectively, by 


(24) 


fb{Eb) ^ ^(Eb-fib)/kBT 


and 


(25) 


ft{Et) = 


1 


\ _|_ IksT ■ 


(26) 


and ksT is the thermal energy. We assume that the Landau levels (LLs) are broadened in energy by F^^t in 
the bottom and top electrodes respectively due to electron - electron interactions, which we model with a set of 
Gaussian functions (to ensure convergence at low magnetic fields) centered on the energies of the LLs (see 
equation]^ [33] 


r(^)= E 




exp 


{E-E^ 

b,t 


(27) 


The density of states is then given by D{E) = (2/7r/^)r(F^). We convert the sum over k states in equation (24) to 
an integral to find the contribution to the current for transitions between LLs rit and is given by 


W{nb,nt) = E_j J \Mbt\‘^dkbdkt, 


(28) 


where L is the device length, so that after using the 6 function to integrate out Et, we find that the current can 
now be expressed by: 


Attc C 

^ = J W{n,,nt) [ME,) - MEt)] D,{E,)Dt{E, - (t>)dE,. (29) 

We model the electrostatics, i.e. the values of and the electrostatic potential energy difference (j), between 
the graphene layers, by solving the following equation: 

(p + Tt) — Pb{pb^ r^) + eV, = 0 (30) 

where d = 1.4 nm is the barrier width, is the charge density on the bottom and top electrodes and the function 
/i(p, F) is found using the density of states, D{E) [5|. From Gauss’s law, and ensuring charge neutrality, we obtain 
the following relationships between V,, V^, p and n,^t' 


e{E,-Eg)=p, (31) 

—eE, = pt, (32) 

where E, = pbjed and Eg = {eVg — p,)/eDg are the fields in the tunnel barrier and gate-oxide barrier respectively 
and Dg = 300 nm is the oxide thickness. 
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Figure 7: Colour maps showing G(Vb, Vg) measured (a) and calculated (b) when B = 2 T and when B = 4 T (c measured, 
d calculated). Colour scales for a,c (b,d) are in fiS (arbitrary units). Filled black circles show loci along which the chemical 
potential in the top and bottom layer, respectively, intersects with the Dirac point in that layer. Lower panels A-D show 
the density of states, D, calculated versus energy, E, in the bottom (red) and top (blue) graphene electrodes and correspond 
to the features labelled A-D in colour maps c and d. Horizontal red and blue dashed lines show position of the chemical 
potentials in the bottom and top electrodes. 


ANALYSIS OF CONDUCTANCE PEAKS 


Fig. 0 shows colour maps of G(V 5 , Vg) = dl/dV^ measured (a,c) and calculated (b, d) when B = 2 T (a,b) and 
B = 4: T (c,d). The parameters used to model the measured data are a = 9 nm and the LL broadening in the 
bottom and top graphene electrodes, F^ and Ft, is set at 4 meV and 4 meV (6 meV and 8 meV) respectively when 
5 = 2 T (4 T). 

In this section we explain in more detail the origin of the conductance peaks observed in G{Vb^ Vg). The filled 
black circles in Fig. [^show the calculated (f4, Vg) loci for which the chemical potential in the top layer intersects 
with the zeroth LL in that layer, see inset A (filled circles running bottom left to top right), and for which 
the chemical potential in the bottom layer coincides with the zeroth LL in the bottom layer, see inset B (filled 
circles running top left to bottom right). Therefore, the local conductance peaks that lie along the X-shaped loci 
correspond to the alignment of the chemical potential in one graphene layer with the peak in the density of states 
for the LL at the Dirac point. 

Fig. shows that in both our experiments and theory, when ^ 5 V and V^ < 0.2 V, increasing 14 initially has 
little effect on G. But when T4 ~ ±0.2 V, there is a sharp increase in conductance. When T4 increases beyond ~ 0.5 
V, G decreases, becoming negative after the peak in /(H). The regions of high G in Fig. [^form stripe patterns 
with similar shapes to the loci marked by the filled circles. This is because they also originate from alignment of 
the chemical potential and LLs when fjib.t = where, in contrast to the yellow curves, rib^t 7^ 0. The crossing 

of these loci gives rise to more islands of high G, for example those labelled “B-D” in Fig. [^c,d. 

When 5 = 4 T, we find good qualitative agreement between the measured and calculated G(14,1^) colour maps. 
Along the loci marked by filled circles in Fig. both maps reveal a series of conductance maxima in similar 
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Figure 8: Colour maps showing comparison of a measured and b modelled G(Vg,Vb) for Vg < 0 and Vb < 0.2 V when 
B = 4 T. Theory curves are calculated with F^ = 3 meV and Ft = 5 meV, cr = 9 nm, and misalignment angle = 1°. 
Circled features and corresponding inset plots show the alignments of LLs in bottom (red) and top (blue) electrodes with 
the chemical potentials indicated by the top of the block colour. 

positions, for example those labelled “B-D” in Fig. and d. As explained above, along the loci, the maxima 
occur as /it sweeps through the LL spectra in the top and bottom layers. The maxima labelled “B” and “C” occur 
when /it coincides with rtt = — 1 and rit = —2 LLs (see insets labelled “B” and “C”). The strength of the maxima 
depends on the alignment of the LLs. For example, the conductance maximum labelled “D” is stronger than “B”, 
because at “B” the LL spectra in the top and bottom layers are aligned and tunnelling occurs from = 0 and — 1 
to Tit = 0 and —1, which have low matrix elements (see main text). By contrast, for case “D” the matrix element 
for tunnelling between the energetically aligned LLs = 3 and = 1 is high. 


Conductance peaks in lower island 


We now analyse the features that appear in the G{Vb^Vg) colour maps at low Vb < ±0.2 V when B = 2 T and 4 
T. These features occur whenever the chemical potential in either the bottom or top layer is aligned energetically 
with one of the LLs in the top or bottom layer respectively. The resulting local maxima in G{Vb^Vg) occur at 
similar positions in the measured (Figs. [^,c) and calculated (Figs. [^,d) colour maps. However, when H = 2 T, 
the theoretical results reveal many more features than the measured data. This is because our calculations assume 
a constant LL width and therefore omit the increased LL broadening that could occur at high Vb in the actual 
device, e.g. due to electron heating. However, the general features of the measured and calculated colour maps are 
similar, in particular the positions of the resonant peaks and the width and shape of the X-shaped low G region. 

In Figs. and b we show an enlargement of Fig. and d focusing on the series of conductance peaks found 
for low Vb and negative Vg when H = 4 T. To model the data at low 14, where electron heating is low, we use a 
narrower broadening (F^ = 3 meV and Ft = 5 meV) than used for the full range of bias voltage. There is very good 
correspondence of the positions of the peaks in the modelled and measured data. As for the local conductance 
peaks considered previously, the peaks arise from a series of alignments of LLs of different index and the alignments 
of the chemical potentials. To aid understanding of these features we highlight two series of resonant peaks labelled, 
respectively, “A-F” and “i-v” showing the alignment of the LLs and the position of the chemical potentials in the 
graphene layers. 
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Figure 9: Colour map showing normalised tunnelling rates, W{nb,nt) (see Eq. 28), for scattering-assisted transitions (taking 
cr = 9 nm) between LLs with indices rih and rit in the bottom and top electrodes calculated using non chiral a and chiral b 
wavefunctions. 


MODEL FOR NON-CHIRAL ELECTRONS 


To understand the effect of pseudospin on our conductance calculations we derive a model for non-chiral electrons. 
The model has the same structure as that presented for chiral electrons but with the electron described by a single 
component wavefunction of the form 




(33) 


where the variables have the same form as those given in section [] Although this form of the wavefunction does not 
correspond to a physical system (it is similar to LL states in III-V materials but with massless Fermions) it allows 
us to distinguish clearly the effect of chirality on the measured and calculated conductance. In Fig. [§L we show 
W{nb^nt) calculated for non-chiral a and chiral b electrons (see Fig. 3c of the main text). The figure reveals that 
for non-chiral electrons, a, transitions between equivalent (c-c and v-v) and different bands (v-c and c-v), have the 
same magnitude; by contrast, for chiral electrons, c-c and v-v transitions are strongly enhanced compared to v-c 
and c-v transitions. 

Fig. [^compares of our conductance calculation for chiral electrons, see section [j with for non-chiral electrons. 
When 5 = 2 T and 4 T within the upper (lower) region, above (below) the dotted and dot-dashed yellow curves, 
the measurements (Figs[TQ^,d) and full calculation (Figs. [T^,f) reveal that the peak amplitudes are largest where 
c-c (v-v) transitions dominate (within the region bounded by the black and white dashed curves). Increasing or 
decreasing outside of this region suppresses the conductance peaks. By contrast, in the calculations using non- 
chiral wavefunctions (Eq. [3^ , the conductance peaks in the lower and upper regions have a constant amplitude 
over the whole range of W (see Fig. [T^,e). This is because the matrix element in the chiral calculations depends 
on the initial and final band of the tunnelling electron and is enhanced for transitions between equivalent bands 
compared to those between different bands. However, in our non-chiral model, the matrix element is equal for 
equivalent transitions between states with the same LL index magnitude for alike and different bands and therefore 
the conduction peak amplitudes are constant across the lower and upper regions. 

A changeover between regions of high and low conductance can also be seen in our recent studies of the G{Vb,Vg) 
characteristics of similar tunnel structures when B = 0. However, in the present work, the changeover is more 
pronounced because the quantizing magnetic field strongly reduces the number of distinct tunnel-coupled states 
that contribute to the current flow; the effect of chirality is strongly magnetic field dependent. Consequently, the 
conductance is more sensitive to the matrix elements for tunnelling between each of these states. 
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Figure 10: Colour maps showing G(Vb,Vg) maps when B = 2 T (a-c) and S = 4 T (d-f) measured a,d and calculated for 
non-chiral (b,e) and chiral (c,f) electrons. Colour scales are in fiS for panels a and d and in arbitrary units for panels b,c,e 
and f. Black and white dashed curves enclose islands within which only conduction-conduction {Vg > 0) or valence-valence 
(Vg < 0) tunnelling occurs. The filled black circles show loci when the chemical potential in the top and bottom layer, 
respectively, intersects with the Dirac point in the corresponding layer. 

Model for non-chiral electrons in zero field 


In zero field we calculate the current for chiral electrons using the model presented in BM- For our calculation 
of current for non-chiral electron, we describe the electrons by plane wave states with the form 

(r) = exp {ikb,t-r) (34) 

where = {kx^ ky) are the wavevectors in the bottom and top electrodes. Therefore the matrix element can be 
found using Eq. 


Mbt = jexp{-a^\kk - kt - AK±|V2). 


(35) 


We then calculate the current by using this form of the matrix element in Eq. (24) summing over fc—states in zero 
field. 
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